###################################################
#
# SECRET WEAPON FOR OIL, ISLAM, and WOMEN PROJECT
#
# Paul Musgrave and Yu-Ming Liou
# musgrave@umass.edu and yl254@georgetown.edu
#
# Last Updated on 27 February 2016
#
##################################################
##################################################
# Section 1: R Housekeeping
##################################################

#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Section 1.1: Clear R workspace
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #

rm(list=ls(all=TRUE))

#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Section 1.2: Set Working Directory 
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #

# Change this to  your working directory
mywd <- "/Users/paulmusgrave/Dropbox/0001 Academic Projects/Ongoing/0127 Oil Islam Women/ISQ Accepted Submission/Replication"

setwd(paste(mywd,"/Data",sep=""))


#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  # 
# Section 1.3: Load Data
#  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #  #
library(foreign)
library(arm)
library(lmtest)
library(sandwich)

longData <- read.dta("ISQFINALRossTSData.dta")


##################################################
# Basic Plot
##################################################

# This is a plot of Michael Ross's model, with the
#  log-squared term


#Clean the Data

secretData <- longData[is.na(longData$FDlaborfemale) == FALSE &
 is.na(longData$FDlogGDP_cap2000_supSQ_1)==FALSE & 
 is.na(longData$FDlogGDP_cap2000_sup_1) == FALSE &
 is.na(longData$FDage_1) == FALSE &
 is.na(longData$FDoil_gas_valuePOPred_1) == FALSE &
 is.na(longData$cty) == FALSE,]

setwd(paste(mywd,"/Drafts/Charts",sep=""))

# Old Version using base
# # RossBase <- function(num){
	# this.num <- subset(secretData,year == num)
	# lm.0 <- lm(FDlaborfemale ~ 
		# FDlogGDP_cap2000_sup_1 + 
		# FDlogGDP_cap2000_supSQ_1 + 
		# FDage_1 + 
		# FDoil_gas_valuePOPred_1,
	# data = this.num)
	# coefs <- summary(lm.0)$coefficients[2:5,1:2]
	# return(coefs)
# }

RossBase <- function(num){
	# Note this version uses coeftest with HC1
	# to duplicate Stata robust standard errors
	this.num <- subset(secretData,year == num)
	lm.0 <- lm(FDlaborfemale ~ 
		FDlogGDP_cap2000_sup_1 + 
		FDlogGDP_cap2000_supSQ_1 + 
		FDage_1 + 
		FDoil_gas_valuePOPred_1,
	data = this.num)
	coefs <- coeftest(lm.0,vcovHC(lm.0,"HC1"))[2:5,1:2]
	return(coefs)
}

years <- sort(unique(secretData$year))
outPutBase <- sapply(years,RossBase)

# Clean the results

incomeCoef <- outPutBase[1,]
incomeSqCoef <- outPutBase[2,]
ageCoef <- outPutBase[3,]
oilRentCoef <- outPutBase[4,]

incomeSD <- outPutBase[5,]
incomeSqSD <- outPutBase[6,]
ageSD <- outPutBase[7,]
oilRentSD <- outPutBase[8,]


# Plot 
par(mfrow=c(4,1),mar=c(3,3,2,1))

# Oil Income Per Capita	
plot(range(years),
     range(0,1.2,-1.2),
     type="n",
     xlab="",
     ylab="Coefficient Value",
     main="Coefficient Estimates of Oil Income Per Capita",
     mgp=c(1.2,.2,0),
     cex.main=1,
     cex.axis=1,
     cex.lab=1,
     tcl=-.1)

abline(0,0,lwd=.5,lty=1)
abline(v=1972.5,lwd=.5,lty=2 )
abline(v=1986.5,lwd=.5,lty=2)

points(years,oilRentCoef,pch=20,cex=.5)
segments(years,
         oilRentCoef-1.96*oilRentSD,
         years,
         oilRentCoef+1.96*oilRentSD,
         lwd=.5)

# we plot Income (normal)
plot(range(years),
	range(0,incomeCoef+1.96*incomeSD,
		incomeCoef-1.96*incomeSD),
	type="n",
	xlab="",
	ylab="Coefficient Value",
	main="Coefficient Estimates of GDP Per Capita",
	mgp=c(1.2,.2,0),
	cex.main=1,
	cex.axis=1,
	cex.lab=1,
	tcl=-.1)
	
	abline(0,0,lwd=.5,lty=1)
	abline(v=1972.5,lwd=.5,lty=2)
	abline(v=1986.5,lwd=.5,lty=2)
	points(years,incomeCoef,pch=20,cex=.5)
	segments(years,
		incomeCoef-1.96*incomeSD,
		years,
		incomeCoef+1.96*incomeSD,
		lwd=.5)

# Now Income Squared
plot(range(years),
	range(0,incomeSqCoef+1.96*incomeSqSD,
		incomeSqCoef-1.96*incomeSqSD),
	type="n",
	xlab="year",
	ylab="Coefficient Value",
	main="Coefficient Estimates of GDP Per Capita Sq",
	mgp=c(1.2,.2,0),
	cex.main=1,
	cex.axis=1,
	cex.lab=1,
	tcl=-.1)
	
	abline(0,0,lwd=.5,lty=1)
	abline(v=1972.5,lwd=.5,lty=2 )
	abline(v=1986.5,lwd=.5,lty=2)
	
	points(years,incomeSqCoef,pch=20,cex=.5)
	segments(years,
		incomeSqCoef-1.96*incomeSqSD,
		years,
		incomeSqCoef+1.96*incomeSqSD,
		lwd=.5)

# Now, a window for age
plot(range(years),
	range(0,ageCoef+1.96*ageSD,ageCoef-1.96*ageSD),
	type="n",
	xlab="",
	ylab="Coefficient",
	main="Working Age",
	mgp=c(1.2,.2,0),
	cex.main=1,
	cex.axis=1,
	cex.lab=1,
	tcl=-.1)

	abline(0,0,lwd=.5,lty=1)
	abline(v=1972.5,lwd=.5,lty=2 )
	abline(v=1986.5,lwd=.5,lty=2)
	points(years,ageCoef,pch=20,cex=.5)
	segments(years,
		ageCoef-1.96*ageSD,
		years,
		ageCoef+1.96*ageSD,
		lwd=.5)


# Save the plot
dev.copy(pdf,'ISQFINALAppendixFigure05.pdf')
dev.off()

##################################################
# Log-GDP Plot
##################################################
# This plot drops the log-squared term

# Reload secretData with new data (no squared term)
secretData <- longData[is.na(longData$FDlaborfemale) == FALSE &
 is.na(longData$FDlogGDP_cap2000_sup_1) == FALSE &
 is.na(longData$FDage_1) == FALSE &
 is.na(longData$FDoil_gas_valuePOPred_1) == FALSE &
 is.na(longData$cty) == FALSE,]



RossLog <- function(num){
	# Note this version uses coeftest with HC1
	# to duplicate Stata robust standard errors
	this.num <- subset(secretData,year == num)
	lm.0 <- lm(FDlaborfemale ~ 
		FDlogGDP_cap2000_sup_1 + 
		FDage_1 + 
		FDoil_gas_valuePOPred_1,
		data = this.num)
	coefs <- coeftest(lm.0,vcovHC(lm.0,"HC1"))[2:4,1:2]
	return(coefs)
}

years <- sort(unique(secretData$year))

outPutLog <- sapply(years,RossLog)

# Clean the results

incomeCoefLog <- outPutLog[1,]
ageCoefLog <- outPutLog[2,]
oilRentCoefLog <- outPutLog[3,]

incomeSDLog <- outPutLog[4,]
ageSDLog <- outPutLog[5,]
oilRentSDLog <- outPutLog[6,]


# Plot -> PDF
par(mfrow=c(3,1),mar=c(3,3,2,1))


plot(range(years),
     range(0,1.2,-1.2),
     type="n",
     xlab="",
     ylab="Coefficient Value",
     main="Coefficient Estimates of Oil Income Per Capita",
     mgp=c(1.2,.2,0),
     cex.main=1,
     cex.axis=1,
     cex.lab=1,
     cex=2,
     tcl=-.1)
abline(0,0,lwd=.5,lty=1)
abline(v=1972.5,lwd=.5,lty=2 )
abline(v=1986.5,lwd=.5,lty=2)
points(years,oilRentCoefLog,pch=20,cex=.5)
segments(years,
         oilRentCoefLog-1.96*oilRentSDLog,
         years,
         oilRentCoefLog+1.96*oilRentSDLog,
         lwd=.5)

plot(range(years),
	range(0,incomeCoefLog+1.96*incomeSDLog,
	incomeCoefLog-1.96*incomeSDLog),
	type="n",
	xlab="year",
	ylab="Coefficient Value",
	main="Coefficient Estimates of GDP Per Capita",
	mgp=c(1.2,.2,0),
	cex.main=1,
	cex.axis=1,
	cex.lab=1,
	cex=2,
	tcl=-.1)
	
	abline(0,0,lwd=.5,lty=1)
	abline(v=1972.5,lwd=.5,lty=2 )
	abline(v=1986.5,lwd=.5,lty=2)
	
	points(years,incomeCoefLog,pch=20,cex=.5)
	segments(years,
		incomeCoefLog-1.96*incomeSDLog,
		years,
		incomeCoefLog+1.96*incomeSDLog,
		lwd=.5)

plot(range(years),
	range(0,ageCoefLog+1.96*ageSDLog,
		ageCoefLog-1.96*ageSDLog),
	type="n",
	xlab="year",
	ylab="Coefficient Value",
	main="Coefficient Estimates of Working Age Share",
	mgp=c(1.2,.2,0),
	cex.main=1,
	cex.axis=1,
	cex.lab=1,
	cex=2,
	tcl=-.1)
	
	abline(0,0,lwd=.5,lty=1)
	abline(v=1972.5,lwd=.5,lty=2 )
	abline(v=1986.5,lwd=.5,lty=2)
	
	points(years,ageCoefLog,pch=20,cex=.5)
	
	segments(years,
		ageCoefLog-1.96*ageSDLog,
		years,
		ageCoefLog+1.96*ageSDLog,
		lwd=.5)

	
setwd("../Drafts/Charts")

# Save the plot
#dev.copy(pdf,'ISQFINALMainBodyFigure2.pdf')


# Hi-res plot
tiff("ISQFINALMainBodyFigure2.tiff", width = 4, height = 4, units = 'in', res = 300)
par(mfrow=c(3,1),mar=c(3,3,2,1))


plot(range(years),
     range(0,1.2,-1.2),
     type="n",
     xlab="",
     ylab="Coefficient Value",
     main="Coefficient Estimates of Oil Income Per Capita",
     mgp=c(1.2,.2,0),
     cex.main=1,
     cex.axis=1,
     cex.lab=1,
     cex=2,
     tcl=-.1)
abline(0,0,lwd=.5,lty=1)
abline(v=1972.5,lwd=.5,lty=2 )
abline(v=1986.5,lwd=.5,lty=2)
points(years,oilRentCoefLog,pch=20,cex=.5)
segments(years,
         oilRentCoefLog-1.96*oilRentSDLog,
         years,
         oilRentCoefLog+1.96*oilRentSDLog,
         lwd=.5)

plot(range(years),
	range(0,incomeCoefLog+1.96*incomeSDLog,
	incomeCoefLog-1.96*incomeSDLog),
	type="n",
	xlab="",
	ylab="Coefficient Value",
	main="Coefficient Estimates of GDP Per Capita",
	mgp=c(1.2,.2,0),
	cex.main=1,
	cex.axis=1,
	cex.lab=1,
	cex=2,
	tcl=-.1)
	
	abline(0,0,lwd=.5,lty=1)
	abline(v=1972.5,lwd=.5,lty=2 )
	abline(v=1986.5,lwd=.5,lty=2)
	
	points(years,incomeCoefLog,pch=20,cex=.5)
	segments(years,
		incomeCoefLog-1.96*incomeSDLog,
		years,
		incomeCoefLog+1.96*incomeSDLog,
		lwd=.5)

plot(range(years),
	range(0,ageCoefLog+1.96*ageSDLog,
		ageCoefLog-1.96*ageSDLog),
	type="n",
	xlab="",
	ylab="Coefficient Value",
	main="Coefficient Estimates of Working Age Share",
	mgp=c(1.2,.2,0),
	cex.main=1,
	cex.axis=1,
	cex.lab=1,
	cex=2,
	tcl=-.1)
	
	abline(0,0,lwd=.5,lty=1)
	abline(v=1972.5,lwd=.5,lty=2 )
	abline(v=1986.5,lwd=.5,lty=2)
	
	points(years,ageCoefLog,pch=20,cex=.5)
	
	segments(years,
		ageCoefLog-1.96*ageSDLog,
		years,
		ageCoefLog+1.96*ageSDLog,
		lwd=.5)

	
	
# Save the plot
dev.off()
	
##################################################
# LogClassical-GDP Plot: Classical Standard Errors
##################################################

# This plot uses classical, not robust, standard
# errors


RossLogClassical <- function(num){
	this.num <- subset(secretData,year == num)
	glm.0 <- lm(FDlaborfemale ~ 
		FDlogGDP_cap2000_sup_1 + 
		FDage_1 + 
		FDoil_gas_valuePOPred_1,
	data = this.num)
	coefs <- summary(glm.0)$coefficients[2:4,1:2]
	return(coefs)
	}

years <- sort(unique(secretData$year))

outPutLogClassical <- sapply(years,RossLogClassical)

# Clean the results

incomeCoefLogClassical <- outPutLogClassical[1,]
ageCoefLogClassical <- outPutLogClassical[2,]
oilRentCoefLogClassical <- outPutLogClassical[3,]

incomeSDLogClassical <- outPutLogClassical[4,]
ageSDLogClassical <- outPutLogClassical[5,]
oilRentSDLogClassical <- outPutLogClassical[6,]


# Plot
par(mfrow=c(3,1),mar=c(3,3,2,1))


plot(range(years),
     range(0,1.2,-1.2),
     type="n",
     xlab="",
     ylab="Coefficient Value",
     main="Coefficient Estimates of Oil Income Per Capita",
     mgp=c(1.2,.2,0),
     cex.main=1,
     cex.axis=1,
     cex.lab=1,
     tcl=-.1)
abline(0,0,lwd=.5,lty=1)
abline(v=1972.5,lwd=.5,lty=2 )
abline(v=1986.5,lwd=.5,lty=2)
points(years,oilRentCoefLogClassical,pch=20,cex=.5)
segments(years,
         oilRentCoefLogClassical-1.96*oilRentSDLogClassical,
         years,
         oilRentCoefLogClassical+1.96*oilRentSDLogClassical,
         lwd=.5)

plot(range(years),
	range(0,incomeCoefLogClassical+1.96*incomeSDLogClassical,
	incomeCoefLogClassical-1.96*incomeSDLogClassical),
	type="n",
	xlab="year",
	ylab="Coefficient Value",
	main="Coefficient Estimates of GDP Per Capita",
	mgp=c(1.2,.2,0),
	cex.main=1,
	cex.axis=1,
	cex.lab=1,
	tcl=-.1)
	abline(0,0,lwd=.5,lty=1)
	abline(v=1972.5,lwd=.5,lty=2 )
	abline(v=1986.5,lwd=.5,lty=2)
	points(years,incomeCoefLogClassical,pch=20,cex=.5)
	segments(years,
		incomeCoefLogClassical-1.96*incomeSDLogClassical,
		years,
		incomeCoefLogClassical+1.96*incomeSDLogClassical,
		lwd=.5)

plot(range(years),
range(0,ageCoefLogClassical+1.96*ageSDLogClassical,
	ageCoefLogClassical-1.96*ageSDLogClassical),
	type="n",
	xlab="year",
ylab="Coefficient Value",
main="Coefficient Estimates of Working Age Share",
mgp=c(1.2,.2,0),
	cex.main=1,
	cex.axis=1,
	cex.lab=1,
	tcl=-.1)
	abline(0,0,lwd=.5,lty=1)
	abline(v=1972.5,lwd=.5,lty=2 )
	abline(v=1986.5,lwd=.5,lty=2)
	points(years,ageCoefLogClassical,pch=20,cex=.5)
	segments(years,
		ageCoefLogClassical-1.96*ageSDLogClassical,
		years,ageCoefLogClassical+1.96*ageSDLogClassical,
		lwd=.5)
	

		
# Save the plot
dev.copy(pdf,'ISQFINALAppendixFigure06.pdf')
dev.off()

